function bin=iflt(img, P, R, nbin)
if nargin < 2
     P = 8;
end
if nargin < 3
     R = 1;
end
if nargin < 4
    nbin = 64;
end
if size(img,3) ~= 1
    img = rgb2gray(img);
end
img = double(img);
[M, N] = size(img);
np = ceil(P/2);
bin = zeros(nbin,4);
% maximum range for mean = -1 .. 1
% maximum range for std = 0..(np/(np-1))^0.5
rng_m = 1.0;
rng_s = (np/(np-1)) ^ 0.5;
u = 0.5*(nbin-1);
pm = u / rng_m; %scale factor
ps = u * 2 / rng_s; %scale factor


[ysize xsize] = size(image);
I = (1:P)* 2*pi/P;
spoints = zeros(P, 2);
spoints(:,1) = -R*sin(I);
spoints(:,2) = R*cos(I);
miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));

% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;

% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));

% Minimum allowed size for the input image depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
  error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end

% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;

% Fill the center pixel matrix C.
C = img(origy:origy+dy,origx:origx+dx);
Ip = zeros(dy+1, dx+1, P);
tic;
for i = 1:P
  y = spoints(i,1)+origy;
  x = spoints(i,2)+origx;
  % Calculate floors, ceils and rounds for the x and y.
  fy = floor(y); cy = ceil(y); ry = round(y);
  fx = floor(x); cx = ceil(x); rx = round(x);
  % Check if interpolation is needed.
  if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
    % Interpolation is not needed, use original datatypes
    N = image(ry:ry+dy,rx:rx+dx);
  else
    % Interpolation needed, use double type images 
    ty = y - fy;
    tx = x - fx;

    % Calculate the interpolation weights.
    w1 = (1 - tx) * (1 - ty);
    w2 =      tx  * (1 - ty);
    w3 = (1 - tx) *      ty ;
    w4 =      tx  *      ty ;
    % Compute interpolated pixel values
    N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...
        w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);
  end  
  Ip(:,:,1) = reshape(N,[dy+1, dx+1, 1]);
end

for i = 1:dy+1
    for j = 1:dx+1
        %Extract I
        I = squeeze(img(i,j) - Ip(i,j,:));
        if ~any(I)
            f = round([u u u u])+1;
        else
            %Normalisation
            I_norm = I / ((sum(I.^2)^0.5)+eps);
            I_norm_1 = I_norm(1:2:end);
            I_norm_2 = I_norm(2:2:end);
            cA1 = (I_norm_1 + I_norm_2) / 2^0.5;
            cD1 = (I_norm_1 - I_norm_2) / 2^0.5;
            %f1->rng1 f2->rng2
            f1 = round([mean(cA1) mean(cD1)]*pm + u); % values in [0..nbin-1]
            f2 = round([std(cA1) std(cD1)]*ps); % values in [0..nbin-1]
            f = [f1 f2]+1;
        end
        s = sub2ind([nbin, 4], f, 1:4);
        bin(s) = bin(s) + 1;
    end
end
toc;
end
